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Abstract 

During the past two decades, our understanding of 
laminar-turbulent transition flow physics has advanced sig- 
nificantly owing to, in a large part, the NASA program 
support such as the National Aerospace Plane (NASP), 
High-speed Civil Transport (HSCT), and Advanced Sub- 
sonic Technology (AST). Experimental, theoretical, as well 
as computational efforts on various issues such as receptiv- 
ity and linear and nonlinear evolution of instability waves 
take part in broadening our knowledge base for this intri- 
cate flow phenomenon. Despite all these advances, tran- 
sition prediction remains a nontrivial task for engineers 
due to the lack of a widely available, robust, and efficient 
prediction tool. The design and development of the LAS- 
TRAC code is aimed at providing one such engineering 
tool that is easy to use and yet capable of dealing with a 
broad range of transition related issues. LASTRAC was 
written from scratch based on the state-of-the-art numerical 
methods for stability analysis and modem software tech- 
nologies. At low fidelity, it allows users to perform linear 
stability analysis and N-factor transition correlation for a 
broad range of flow regimes and configurations by using 
either the linear stability theory (LST) or linear parabolized 
stability equations (LPSE) method. At high fidelity, users 
may use nonlinear PSE to track finite-amplitude distur- 
bances until the skin friction rise. Coupled with the built-in 
receptivity model that is currently under development, the 
nonlinear PSE method offers a synergistic approach to pre- 
dict transition onset for a given disturbance environment 
based on first principles. This paper describes the govern- 
ing equations, numerical methods, code development, and 
case studies for the current release of LASTRAC. Practi- 
cal applications of LASTRAC are demonstrated for linear 
stability calculations, N-factor transition correlation, non- 
linear breakdown simulations, and controls of stationary 
crossflow instability in supersonic swept wing boundary 
layers. 
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Introduction 

In the late 80's till the early 90's, laminar-turbulent 
transition research was one of the outstanding issues for 
the National Aerospace Plane (NASP), High Speed Re- 
search (HSR), and Advanced Subsonic Technology (AST) 
programs. After a dormant period, there has been re- 
newed interest in supersonic laminar flow aircraft in recent 
years due to the DARPA Quiet Supersonic Platform (QSP) 
project. The main objective of the QSP project is to eval- 
uate the feasibility of a low-boom laminar flow supersonic 
aircraft.^ The ongoing NASA Langley Supersonic Vehi- 
cle Technology (S VT) program is also aimed at developing 
technologies required for designing such supersonic vehi- 
cles. Computational tool development as well as a series of 
wind-tunnel and flight experiments are planned. 

Laminar flow control has been investigated extensively 
over the past decades. Experimental, theoretical, as well 
as computational efforts on various issues such as receptiv- 
ity, linear and nonlinear evolution of instability waves take 
part in broadening our knowledge base of this intricate flow 
phenomenon. For low speed flows, various stages of the 
transition process are now fairly well understood. Theory 
and numerical simulations using either parabolized stabil- 
ity equations (PSE) or direct numerical simulations (DNS) 
accurately describe disturbance generation and evolution 
until the transitional stage. Secondary instability in two- 
dimensional and swept wing boundary layers observed in 
low-speed experiments may also be calculated accurately 
using nonlinear PSE and the eigenvalue approach.^' ^ Saric 
ei al."^ discovered that stationary crossflow instability can 
be controlled by introducing distributed roughness near the 
leading edge to alter the mean state and subsequently sup- 
press or delay the growth of the most unstable disturbances. 
It has been shown that this control phenomenon may be cal- 
culated and parametrically studied by nonlinear PSE (see 
Malik et al.^ and Chang'^). 

For supersonic to hypersonic boundary layers, the break- 
down mechanism is not fully understood largely due to the 
lack of experimental measurements. For two-dimensional 
boundary layers, it has been shown computationally that 
the secondary instability mechanism observed in low speed 
cases may still be present.^ But other mechanisms, such 
as oblique-mode breakdown (Chang & Malik^), may be 
more relevant because of the large growth of the oblique 
disturbances in low supersonic boundary layers. At hy- 
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personic speed, large amplitude second-mode disturbances 
have been observed experimentally^'^ but the breakdown 
mechanism is unclear, mainly because the wave orienta- 
tion is largely unknown. The recent experiment of a flared 
cone in a conventional Mach 6 tunnel by Horvath et al.^^ 
indicates that transition may or may not be caused by 
second-mode disturbances because of a small N-f actor at 
the transition onset. We need more detailed experimen- 
tal measurements to clarify this issue. Receptivity of hy- 
personic boundary layers to acoustic disturbance has been 
studied by Federov^^ and Zhong.^^ The role of the entropy 
layer on hypersonic boundary layers is investigated both 
theoretically and numerically by Fedorov & Tumin.^^ For 
supersonic swept wing boundary layers, Saric is currently 
extending his distributed roughness control for supersonic 
flows. 

Despite all the advances in transition research, the tra- 
ditional way of transition correlation and prediction using 
the f method remains one of the mainstream tools used in 
the design phase due to its simplicity and availability. The 
N- factor method provides a simple way to correlate transi- 
tion locations and in the meantime predict transition onset 
using a prescribed N-value. For a higher fidelity transi- 
tion prediction based on first principles, one can start from 
the receptivity model and compute the disturbance initial 
amplitudes according to the given disturbance environment 
such as a roughness or free-stream disturbance distribution. 
Receptivity calculations may be based on a simple lin- 
ear theory or a more advanced linear Navier-Stokes (LNS) 
method. The evolution of the internalized disturbances 

with known amplitudes can then be calculated using non- 
linear PSE all the way until the skin friction rise. Alterna- 
tively, a secondary instability calculation may be performed 
for a steady, highly nonlinear boundary layer with the pres- 
ence of stationary crossflow vortices.^ To further refine the 
transition region or to carry the perturbed field into a fully 
turbulent stage, DNS or large eddy simulation (LES) may 
be used. 

The integrated transition prediction approach described 
above encompasses several mature computational means to 
form a package of tools for high-fidelity predictions beyond 
the conventional method. Our current research effort 
is geared toward this integrated approach. In addition to 
the tool development, at issue is the disturbance environ- 
ment. A complete description of the environment remains 
very difficult except for a forced transition where manu- 
ally excited roughness or free-stream waves outweigh the 
natural ones and thus may be quantified accurately in the 
computations. As a result, statistical means must be intro- 
duced to account for the natural environment. Some efforts 
have been made along this direction.*^ A case study for su- 
personic swept wing using the integrated approach is also 
available. 

One of the main goals of the ongoing NASA Langley 
projects for transition flow physics and prediction is to de- 
vise a set of tools to enable traditional as well as integrated 
transition prediction. The Langley Stability and Transi- 


tion Analysis Code (LASTRAC) was developed to meet 
two main objectives: to provide an easy-to-use engineering 
tool for routine use and to incorporate state-of-the-art com- 
putational and theoretical findings for integrated transition 
predictions. For the former objective, LASTRAC can per- 
form linear calculations and transition correlation using the 
N -factor method based on either LST or a more advanced 
linear PSE method. The user interface of LASTRAC code 
was designed to allow users to identify unstable parameter 
space in terms of disturbance frequency and wave numbers 
with minimum effort. Transition prediction or correlation 
thus can be made with the identified parameter space. In 
addition to the traditional LST-based N-factor calculations, 
LASTRAC allows users to perform linear PSE N-factor 
calculations, which in some cases give a more compact N- 
factor range and better correlation. 

To achieve the second objective, LASTRAC provides 
transition simulation capability based on an absolute am- 
plitude. The receptivity module computes the initial am- 
plitudes near the neutral stability location. Nonlinear PSE 
then simulates disturbance evolution until skin friction rise. 
For further refinement near the transition region, a sub-grid 
scale model based nonlinear PSE, or large eddy simulation 
(LES) and direct numerical simulation (DNS) may be used. 
These modules are all currently under development. We 
will discuss mainly the nonlinear PSE method in this pa- 
per. 

In the current LASTRAC release, LST and linear and 
nonlinear PSE options are implemented for 2-D, axi sym- 
metric, and infinite swept wing boundary layers. In addi- 
tion, 2-D mixing layer, axi symmetric jet and vortex flows 
are also handled by the current release. Extension to gen- 
eral three-dimensional(3-D) boundary layers is currently 
under development. This paper is a status report of the 
LASTRAC development. We discuss governing equations, 
numerical solutions, code development, and case studies of 
LASTRAC for practical examples in this paper. 

Governing Equations 

The problem of interest is a compressible flow over a 2- 
D or axi symmetric body with a coordinate system depicted 
in Fig.l where x is the streamwise, y the wall-normal, 
and c the span wise direction (or the azimuthal direction 
for axisymmetric cases). For the case of infinite swept 
wing boundary layers, the coordinate system used is shown 
in Fig. 2 where x is the normal-chord direction and c is 
parallel to the leading edge. A nonorthogonal coordinate 
system can also be used for infinite swept wing boundary 
layers and will be used for future LASTRAC releases. A 
body -fitted coordinate is adopted because we are interested 
in a region adjacent to the wall. For 2-D or axisymmet- 
ric configurations, both streamwise curvature along x and 
transverse curvature along c may be included in the compu- 
tation; only streamwise curvature is allowed for an infinite 
swept wing configuration. In the orthogonal body-fitted 
coordinate system, elements of length are h\dx, dy, and 
hsdz, in the x-, y-, and c-direction, respectively. The 
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element length can be written as 

(II = \/{h\(h')~ + (lir + 

The metric h i is associated with the streamwise curvature 
K and is defined as 

h \ ~ 1 + 

The metric h-^ is unity for nonaxi symmetric bodies; for ax- 
i symmetric configurations, it is defined as 

A 3 = Vh + yCO>'-<{0) 

where is the local radius and 0 is local half-angle along 
the axi symmetric surface. 


Shock 



Fig. 1 Coordinate system for a 2-D or axisymmetric bound- 
ary layer 

The evolution of disturbances is governed by the com- 
pressible Navier-Stokes equations 
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Fig. 2 Coordinate system for an infinite swept wing boundary 
layer 

The perfect gas relation, 

p ^ p'^T. (4) 

is used as the equation of state. 

Equation 1 is in a nondimensional form, i.e., all lengths 
arc scaled by a reference length scale ( , velocity by ih , den- 
sity by p^ , pressure by pt ui, temperature by 7^ , viscosity 
by /le. and time by l/ii , . For stability calculations, these 
normalization quantities are taken to be the boundary-layer 
edge values. The length scale ( is the similarity boundary- 
layer length scale defined by 

( = 

and the corresponding Reynolds number is defined by 

R = a, i/ue . (6) 

The solution to the Navier-Stokes equations consists of 
two parts, the mean laminar flow solution and the distur- 
bance fluctuation, i.e.. 


where V is the velocity vector, p the density, p the pres- 
sure, 7 the temperature, specific heat, k the thermal 

conductivity, and p and A the first and second coefficient of 
viscosity, respectively. The viscous dissipation function is 
defined as 


d> = A(v ■ V')- + |[vr + vr^]^ (2) 

The second viscosity is related to the first one by the Stokes 
parameter s, defined as 


u — u -f i/. r = V -f I ^ u' = U’ -f iv^ 

p — jt ^ p\ p — p p' . T — T -h T' (7) 

p — -h /y^ A = A + A^ A' = A* -h 

Substituting Eq.7 into Eq. 1 and subtracting the governing 
equations for the steady-state basic flow, we obtain the gov- 
eming equations for the disturbances as 


do ^4 do ^do c do 
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I ,Vr.r0-O Trjy (‘^“O .. (}~0 

u ( ~T^ — f"^ Vf/TpT (8) 

no n~\ (h-~ hi ().r(}y ' ' (h/~ 

\\ r : O ~0 \ ()-0 \ , ()'0 

h:i (),v()z hii dijc): //rj ():- ^ 

where 0 is the disturbance vector defined as 

o = (;/. u\ v', irW)^. (9) 

The coefficient matrices, l\ A, B, f\ I), ^!ry» etc. 

are the Jacobians of the flux vectors, and /?o = Ut (\)/iy^ is 
the Reynolds number used to normalize the equations. The 
reference length scale fo depends on the flow configura- 
tion. For a boundary-layer flow, the boundary-layer length 
scale (defined in Eq.5) computed at the initial location of 
the stability calculations is used. For a jet or shear layer, 
the jet diameter or shear layer thickness may be used. 

The Jacobian matrices can be further decomposed into 
two parts. The first part contains only mean quantities and 
the second part contains perturbation quantities. For exam- 
ple. 


A = .4^ -h .4^^ 

where .4^ = A^(p. r/. r, f, ...) and .4” 

For a 2-D or infinite swept boundary layer, we assume 
the disturbance field is periodic in both temporal and span- 
wise directions. The disturbance vector 0 can then be 
expressed as the following discrete Fourier series: 

M N 

0{x.y.zj)= ^ (10) 

m = ~M n = -N 


mode shape of a mode (h?. n) is represented by the com- 
plex vector \ 

Equation 8 is the governing equation for the disturbance 
evolution. Direct solution of Eq.8 is referred to as the di- 
rect numerical simulation (DNS) method. DNS requires a 
significant amount of computational time even for a simple 
linear case. Hence, a more efficient approximate solution 
to Eq.8 is more desirable. When the disturbance amplitude 
is very small relative to the mean quantities, each Fourier 
mode evolves independently. Nonlinear interaction among 
disturbance modes is negligible. The linearized version of 
Eq,8 can be written as 


1 00 A^ O 0 ,00 (0^ 00 , 


Oi/ h:i Oz 


1 /^TrO~0 , Kry O~0 _ , . O~0 


(-^—4- -f 

D ' 1.2 ' 


Hi) h\ Ox- hi Oxiiy 


rr + y, 


' (hr 


(14) 


vj.. d-d> Vy'' d-0 0-0 

I 13 dxdz /?;j dydz hi dz- ' 

For a single disturbance mode (u.\ /:(), 0 can be expressed 
as 

0 = \(j-,.v)f’ (15) 
Substitute the above equation into Eq.l4, we obtain 


,-■1' mL.dx 


)ir + Uy + 




h:iRi, ()x 

, „/ . ...I ii3C ^ 

+(D — /w'l 4 — 1- 


h^Ri) dy 
ifV!. ^ 

)\ = (16) 


h'iRo 


R,r h-i dx~ hi dxdy 


+ 1 


■' — ) 


where lj and A are the fundamental temporal and span- 
wise wave number of interest, respectively. M and N are 
the number of Fourier modes included in the calculation 
in time and space. The wave number w is related to the 
physical frequency / by 


Another nondimensional frequency often used in stability 
calculations is defined by 


Eq.l6 is the linearized Navier-Stokes (LNS) equation, 
which is a set of constant coefficient PDF’s. Solution of the 
LNS still requires a significant amount of computational 
time. To further reduce it, approximations to Eq.l6 must 
be made. 

In the PSE approach, as suggested by Herbert and 
Chang et al.,^ we decompose the mode shape \{ j\ y) into 
two parts; the wave part described by a complex wave num- 
ber a and the shapefunction part V' to track the nonparallel 
effect, i.e., 


F = 


27ny, 

, ,0 / • 


The span wise wave number is defined by 


(12) 


0=2t:/X, (13) 

where is the spanwise wavelength normalized by the 
length scale 4. The fundamental mode (u;, 0) determines 
the size of the computational domain( 1 // in time and A^ 
in the spanwise direction) and A/ and A' represent the nu- 
merical resolution (2A7 -f 1 and 2 V -i- 1 discretized points 
in time and the spanwise direction) in the simulation. The 


\(j',//) = il’(x.y)e (17) 

in which 0 = (;>.?/, r. u , T). Here, the integral form 
is used for the wave part to allow a variation in x and 
to record the history effect of o. The shape function V 
varies both in x and y. Hence, nonparallel effect is re- 
tained through the variation of stream wise wave number 
and shape function. Substituting Eq, 17 into Eq.l6, we have 



+B^+D0 

ay 


. yly O' t , 

Ri, h'i dx^ hi dxdy'^ ^«dy'^ 

(18) 
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where matrices .1, B, and 1) are defined by 


.1 

B 

D 


.1' 2m in’ 


-h 


1 

h'^Ri) 

} /o 1 rf/ 



llalii) 


inA^ i> 

— iuj r' -|- 

/n ^ 



■ ■ 4- 

— r . 


do 


(19) 


hih-Aiia /^r;/?n 


= - 
+ 
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f)/ h\ c).r (iy If-.i (F 


1 

7 ^ 


y/j r)j'- 


_L \ 

//i dj‘()y ()y- 




+ 


wwy-c^ 


/^ I /?3 (hF^' 


+ -rt-77^) 


/)r^ ryc- 


where the disturbance vector o is defined in physical space 
(see Eq,9). Nonlinear forcing for each Fourier mode is de- 
fined by the following Fourier transform: 


The above PDF differs from the LNS equation (Eq. 1 6) only 
in the terms associated with the wave part o . In fact, a wave 
part may also be introduced in the LNS solver to increase 
the numerical accuracy. The harmonic LNS solution pro- 
cedure may be used to solve Eq.18. However, to improve 
computational efficiency, a marching procedure is more de- 
sirable as long as Eq.18 is parabolic in .r. To this end, we 
neglect the viscous derivatives in x and Eq.18 is reduced to 


.4 


dv' 

dx 


+ + Or = 

(>y 


Hi) (>y- 


( 20 ) 


Chang et al.- pointed out that Eq.20 is only “nearly 
parabolic” due to the streamwise pressure gradient term 
0pldi ) inherited from the original Navier-Stokes equa- 
tion 1. Numerical instability similar to that observed in 
the parabolized Navier-Stokes (PNS) equations'^ will oc- 
cur if one attempts to solve Eq.20 by using a small enough 
marching step size. They further suggested that Vigneron's 
approximation'^ may be used to suppress the numerical in- 
stability. Li and Malik’" use Fourier analysis to prove the 
existence of numerical instability and quantify the numeri- 
cal instability bound. 

If we retain all nonlinear terms, the governing equations 
become 


rv = -M = 

(24) 

The forcing term Fmn for a given Fourier mode u ) can 
be evaluated by collecting all nonlinear terms contributing 
to it. For instance, the (0. 1) mode interacting quadrati- 
cally with the (1.0) mode would contribute to the (1,1) 
mode. In LASTRAC, nonlinear forcing terms are evaluated 
by computing in physical space and then transforming 
F’ ’ back into the wave number space to obtain • 

The governing linear and nonlinear PSE are solved by 
a tburth-order central difference in the wall -normal direc- 
tion and a first- or second-order one-sided difference in 
the marching direction. By invoking the locally parallel 
assumption, the linear PSE equations may be reduced to 
the LST eigenvalue problem. LASTRAC is capable of 
handling LST and linear or nonlinear PSE solutions for 
a user-supplied mean flow. Details of the discretization 
scheme and boundary conditions may be found in Ref. 
The PSE marching solution requires an initial condition, 
which is discussed in the next section. 

Initialization of PSE Marching 


, dl'mn , D . n f 


( 21 ) 


where matrices , , and are defined by 


-477} ri — I 


. 4 ^ 2/0 77)7) in 3 \\ 


hi hfRo h:^Ro 

„/ 1 J-y _ iuJ\ p. 

, j /o77,7).4^ indC^ 

D„„, = D‘ - wiuir' -I- — 7 -I- — r 


( 22 ) 


^ no3V^ ^ 


dx 


hifisRo hr^Ru 


and .477)7) = The left-hand side of Eq.21 con- 

tains only linear coefficient matrices. All nonlinear terms 
are included in the forcing function Fmn^ which is the 
Fourier component of the total forcing, F" , defined as 


The PSE formulation requires an initial condition to start 
the streamwise marching procedure. A simple way to ini- 
tiate marching is to use the eigenfunction obtained from 
the quasi-parallel LST solution. However, because of the 
neglected nonparallel terms, the marching PSE solution 
obtained this way would have a “transient" region where 
the PSE solution attempts to slowly adjust the inconsis- 
tent parallel eigenfunctions. This transient effect, more 
noticeable for supersonic boundary layers, is a result of the 
difference in governing equations between LST and PSE. 
To alleviate the transient effect, one can use the multiple- 
scale method.^"’ The computed nonparallel eigenvalue 
and eigenfunction may then be used to initiate the PSE 
marching procedure. 

LASTRAC offers an alternative nonparallel eigen- 
solution (NES) procedure based on the discretized parab- 
olized governing equations. The governing equation is 
transformed to a general coordinate system by using 
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Numerical calculation is then performed in the (^ . //) plane. 
Using first‘Order finite difference in ^ , the governing linear 
PSE at two consecutive solution locations, i and / + I , may 
be written as 


-h — l’/) -h V 

()lj 

D,+ iv-,+ ^ + .i, + i(r,+ i - (,',)+ S, + i^%ti(25) 

ai] 






rhf^ • 


For this equation, we assume that the second derivative in 
^ is negligible and Oi /Oi, is the same for both point i and 
/ + 1. This assumption is equivalent to taking forward dif- 
ference for point i and backward difference for point / -f L 
The stream wise wave number a, and is related by 


da 

Ojr' -j- — 

dx 


(26) 


in which Ax is the stream wise distance between points i 
and 7 -f 1 . 

If we define the combined shape function vector ^ as 




V'i 


then Eq.25 may be rewritten as 


-dMf 


where matrices D, B, and V are 


(27) 


D 

B 

V 


Di — Ai A I 

4 , 1 77/ ^ j -4 j _j_ ] 

Br 0 \ 

0 ^.+ 1 ; 

^ 0 4 

0 + l J 


(28) 


At the boundaries, homogeneous boundary conditions (e.g. 
no slip at wall and Dirichlet or nonreflecting conditions at 
free stream) may be applied for both ? and /+ 1 stations. Eq. 
27 is a set of homogeneous ordinary differential equations 
(ODE). After discretization in t], Eq. 27 can be cast into 
the following nonparallel eigenvalue problem : 


= 0 (29) 

where = Lnp(<^.a,da/dxjA, For a given distur- 
bance (wu', /:^), two eigenvalues a and da/dx need to be 
found. Unlike the LST eigenvalue problem, the NES has 
one more unknown: da/dx. Since the variation of a is 
generally small, our first approximation is that da / dx = 0 . 


Thus the eigenvalue problem is similar to the LST case, ex- 
cept that the operator now contains solutions at two 
consecutive streamwise stations and requires more time 
to solve (for both local and global eigenvalue searches). 
In fact, obtaining the global spectrum with both n and 
da/d.v as unknowns is nontrivial. In LASTRAC, we as- 
sume f/n /dx = 0 for the solution of the global eigenvalue 
spectrum. 

For the local nonparallel eigenvalue search, as men- 
tioned above, we may assume o, = o ,+i. The solution 
procedure is then similar to the LST case. Alternatively, 
we may use a local Newton's iteration to solve for n , and 
0 / 4.1 simultaneously. To this end, we need two target con- 
ditions to be iterated upon. If we use u-velocity at the wall, 
then the wall boundary condition i/j(U) — = 0 

is replaced by p/( 0 ) = p, 4 . 1 (G) = 1 . This implies that 
we will normalize the eigenfunctions such that the pressure 
eigenfunction is unity at the wall. The relaxed u-velocity 
condition at the wall becomes the target of the Newton's 
iterations. In this case, we further assume that the pressure 
shapefunction at the wall does not vary from station i to 
7 + 1. For most cases, the wall pressure variation dp/dx is 
small and the above approximation is valid. If /] = 7 //( 0 ) 
and /'J = 4. 1 ( 0 ) , using first-order expansion, we have 


G — U 4- -h ( ^ (30) 

0 = h + + l 

uaj 

where Ao, = q"+' - of and Ao, + ] = q"+/ - q'’^, 

denote the change of eigenvalues from iteration n to itera- 
tion n -j- 1 . Solution of Eq. 30 provides two eigenvalues: 
and 0 ; 4 _]. The eigenvalue iteration is repeated until 
Ao, and Aq , 4 _i are both smaller than a prescribed toler- 
ance. We have also used the wall temperature instead of the 
u-velocity as the target condition for Newton's iterations; 
numerical experiments indicate that both targets worked 
well. For hypersonic boundary layers, relaxing wall tem- 
perature tends to be more robust than wall velocity. 

The above non-parallel eigenvalue formulation has sev- 
eral distinct characteristics. First, for every eigensolu- 
tion, we obtain eigenfunctions at two consecutive loca- 
tions. Therefore, streamwise derivatives and thus the ef- 
fective nonparallel growth rates accounting for this stream- 
wise variation may be evaluated. Second, local eigenvalue 
search using only one alpha (da/dx = 0 ) is more robust 
because the global eigenvalue solver is a good approxi- 
mation to the local solver. However, due to the neglected 
da/dx, this approximation produces a small but noticeable 
transient effect (smaller than using the LST solutions) if we 
use the obtained eigenmode to initiate the PSE marching 
procedure. Again, this transient effect is more pronounced 
in supersonic than in low-speed boundary layers for 2-D 
flows. For crossflow instability, the variation of a near 
the leading edge is usually large and the one-alpha proce- 
dure yields a noticeable transient effect even at low speed. 
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Fig. 3 Linear PSE solution obtained by using three different 
initial eigenfunctions, compared with the solution initiated far 
upstream to assess the transient effect 

Furthermore, even though this transient effect is noticeable 
in the disturbance growth rate, it is much less noticeable 
for the disturbance amplitude which represents the inte- 
grated effect of the growth rate. Figure 3 compares linear 
PSE solutions obtained by the one-alpha, two-alpha, and 
quasi-parallel LST eigenfunctions. The solution initiated 
far upstream is also plotted as a baseline solution. Clearly, 
the two-alpha eigenfunction produces the least transient ef- 
fect. 

On the other hand, the two-alpha iteration procedure is 
less robust because it requires convergence of two eigen- 
values and only first-order Newton’s method is employed. 
Deriving a second-order method is possible but would re- 
quire introducing more target conditions. The lack of ro- 
bustness also stems from inconsistency between global and 
local eigenvalue solvers (noted that in the global solver, 
we assume da/dx — 0). The initial guess of do/dx is 
important for the two-alpha procedure. We have found 
that for 2-D boundary layers, using the guess from the 
global solver and imposing da/dx = 0 generally leads to 
a converged solution. However, for crossflow instability 
near the leading edge of a swept wing, providing a good 
guess of da /dx is crucial to convergence of the eigenvalue 
search. LASTRAC uses several built-in da/dx guesses for 
the two- alpha eigenvalue search. However, in some cases, 
users need to provide a good guess in order to get a con- 
verged eigenvalue. It was also found that the eigensolution 
obtained with the two-alpha procedure produces minimal 
transient effects when used to initiate the PSE marching. 

As discussed in the previous section, the determination 
of is in fact nonunique within the nonparallel frame- 
work. The above nonparallel eigenvalue formulation typ- 
ically generates more discrete eigenmodes than its quasi- 
parallel LST counterpart. Therefore, selecting eigenmodes 


of interest from the global eigenvalue spectrum becomes 
more difficult. Users may need to provide certain addi- 
tional criteria to filter out spurious modes. Interestingly, 
one quasi-parallel eigenmode would split into two modes in 
the corresponding nonparallel global spectrum mainly be- 
cause we are solving two consecutive eigensolutions simul- 
taneously and we assume that a, = 0 ;+i. We found that 
with the two-alpha local eigenvalue search, two split modes 
converge to the same mode. Mathematically, the nonpar- 
allel eigenvalue formulation introduces more degrees of 
freedom in the system and thus allows more discrete eigen- 
modes (other than the continuous spectrum). The mode 
closely associated with the unstable mode given by the par- 
allel LST still appears as an eigenmode in most cases. Near 
the leading edge of the boundary layer, instability wave 
formation and evolution is governed by the receptivity pro- 
cess; therefore the nonparallel eigenmode may or may not 
be relevant from the theoretical standpoint. LASTRAC se- 
lects the most unstable mode that gives the largest growth 
rate at the location of interest to initiate the marching pro- 
cedure. Numerical exp>eriments indicate that the eigenvalue 
and eigenfunction generated by the NES (both one- and 
two-alpha approaches) produce a smaller transient effect 
than those from the quasi -parallel LST solution. 

Code Development 

The LASTRAC code was developed from scratch. Sev- 
eral guiding principles were set forth in the initial stage of 
L/\STRAC development. They are: 

• Object-oriented (OO) design and implementation 

• Generic programming using templates 

• Optimization to avoid abstraction penalties and ensure 
efficiency 

• Multithread and message passing interface (MPI) 
based parallelization 

• Source code control and some configuration manage- 
ment 

The first two items limit the choice of programming lan- 
guage to C-HH- because it is the only mature language that 
supports both object-oriented and generic programming. 
The OO paradigm has been recognized by the software 

orld as a mainstream approach for software development. 
An OO software system may take longer to develop but it 

ill be much easier to maintain and extend as the project 
goes along. The concepts of data abstraction, inheritance, 
and polymorphism are the themes of an object-oriented 
language. Separation of interface and implementation prin- 
ciple makes software modules more independent of one 
another. 

The common concept for modem software development 
is to divide the software into a class of similar but not iden- 
tical objects with a common interface. The clients of these 
objects only need to know about the interface — not the 
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implementation details. In this way, new objects may be 
derived and added to the system without needing to change 
the common interface or any existing object client's code. 
Code reuse is thus achieved via reciting common interfaces 
for various class objects but not through cut-and-paste of 
existing codes. 

The OO design procedure typically employs a formal 
requirement analysis along with use case studies. The Uni- 
fied Modeling Language (UML)^*^ is a widely used tool for 
OO analysis and design. Due to the lack of resources, we 
did not pursue the UML route for analysis and design. In- 
stead, the structure of LASTRAC was laid out based on 
the numerical formulation and then different design pat- 
terns were applied when appropriate. The design patterns 
used in LASTRAC came mostly from Gamma et al.^*’ with 
some modifications to cope with numerical efficiency. A 
design pattern is a set of objects and interfaces that occur 
over and over again in a software system and hence is also a 
solution to a software problem in a context. Many patterns 
proposed by Gamma et al. are used in LASTRAC. For ex- 
ample, the Singleton pattern is used for input and control 
parameters, a Factory Method is used to generate objects 
on the fly, the Observer pattern is used to let the metric 
class, coordinate transformation class, and Jacobian matri- 
ces class update their contents when the solution marches 
to the next station. The Template Method and Command 
patterns are used to implement a multithread class that al- 
lows a nonstatic class member function to be called within 
a C-style thread routine. 

Generic programming^^ is a new paradigm that allows 
a programmer to introduce a parameterized class mem- 
ber or member function in a function or class definition. 
The introduced parameter is called a template. A routine 
designed this way allows unrelated classes or algorithms 
to be plugged in as a template parameter as long as the 
plugged-in classes fulfill the compilation requirement. For 
example, one can design a numerical integration class that 
allows various approximation rules be specified as a tem- 
plate parameter. The C+-»- programming language includes 
a powerful library called standard template library (STL) 
based on the generic programming concept. LASTRAC 
uses templates extensively in the class design. STL con- 
tainers and algorithms are used to ensure efficiency and 
flexibility. 

Discussing the details of the LASTRAC code design and 
structure is beyond the scope of this paper. Here, we only 
describe important elements of the code structure. Users 
of LASTRAC need to prepare two files: one contains the 
mean flow; the other has the input parameters instructing 
which computation (LST, linear or nonlinear PSE, etc.) 
needs to be performed. The input parameter file is read 
by LASTRAC through a parser. The parser takes an in- 
put format similar to the namelist format in Fortran but 
allows C++ style comments. LASTRAC performs a “san- 
ity check" for the input parameter file and flags errors if any 
user mistake is found. The mean flow reader was designed 
to be a hierarchy of classes that can be extended for per- 


fect gas/reacting flow, 2-D/3-D through inheritance of the 
abstract class. 

The core part of the code uses the observer pattern (see 
Gamma et al.^‘') where two subjects, LocalSubject and 
MarchingSubject, perform a local eigenvalue solution and 
a marching PSE solution, respectively. The observers are 
the coordinate transformation, metric and derivatives, and 
Jacobian classes. These classes update their states when the 
LocalSubject and MarchingSubject finish their calculations 
and move on to the next location. The nonlinear march- 
ing subject communicates with LocalSubject and March- 
ingSubject to obtain a local eigenvalue for initiating the 
marching and to solve for each individual mode during the 
nonlinear iteration. It has several logistic classes, such as 
a Fourier container class, to perform fast Fourier transform 
and to handle Fourier-mode related operations. 

The finite difference operation and boundary condition 
treatment is handled by using a template class and the 
generic programming paradigm. These generic classes 
would work for different orders of differentiation and vari- 
able types such as real or complex. The left hand side 
operator of the PSE and eigenvalue solver were cast in an 
operator form to facilitate future maintenance. 

To achieve computational efficiency in the LASTRAC 
code design, we tried to avoid using many fine-grained 
objects. The use of templates in generic programming al- 
lows substitution of template parameters at compile time; 
consequently the code can be optimized for run-time per- 
formance. In LASTRAC, the template meta-programming 
technique is used to optimize the matrix multiplication pro- 
cedure and the block matrix solvers. Extensive use of 
generic programming and template utilities makes LAS- 
TRAC as efficient as an earlier code written in Fortran 77. 

Parallelization for the nonlinear PSE option in LAS- 
TRAC is implemented by using the multithread technique 
in conjunction with the more traditional MPI approach. 
Posix thread is used in LASTRAC to guarantee portabil- 
ity. Multithread and MPI options are only implemented 
for nonlinear PSE calculations. Linear options, including 
LST and linear PSE, can only be run as single-thread ap- 
plications. Since linear calculations are mode independent, 
users can launch many LASTRAC runs concurrently on a 
multiprocessor or clustered environment. 

For nonlinear PSE calculations, each Fourier mode is 
solved by a thread or a process under MPI. However, the 
nonlinear forcing terms have to be computed collectively. 
Parallelization of nonlinear forcing is only done through 
computing terms in physical space simultaneously. Each 
MPI process may spawn multiple threads. For a multi- 
processor node in a cluster, this combined thread and MPI 
approach makes LASTRAC run more efficiently. 

Source code control and configuration management is 
performed through the public domain software CVS. Unit 
testing was accomplished for major modules during the ini- 
tial development stage. Due to lack of resources, we do not 
follow unit and integration testing procedures rigorously. 

A semi-automated integration testing procedure is also im- 
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plemented. 

LASTRAC code has been ported to all major Unix plat- 
forms. In the release, we include executables for Linux 
Intel, SUN Solaris, and SGI. We are planning to port the 
code for Linux Alpha and Windows in the near future. 

Code Validation and Test Caj^ 

The LASTRAC code was validated by comparing its re- 
sults with those obtained in the literature and by other pub- 
licly available codes. The first case used for validation was 
the quasi-parallel LST results published by Malik." The 
test cases range from low-speed (incompressible) to hyper- 
sonic conditions. Linear and nonlinear PSE results were 
compared with those given in Chang et al." "^ "^ For swept 
wing boundary layers, the Arizona State University test 
configuration^^ was used for validation, and results were 
compared with Malik et al.^ For all these validation cases, 
the results compared very well with existing solutions. To 
demonstrate the use of LASTRAC, we discuss several test 
cases in details below. 



Fig. 4 Maximized N-factor vs. Reynolds number for / from 
2 kHz to 45 kHz for a Mach 2 flat plate boundary layer 

Mach 2 Flat Plate 

The first case is a Mach 2 flow over a flat plate. LAS- 
TRAC reads the mean flow and prints out a summary as 
follows: 

/* Mean Flow Parameters Reading in *\ 



Fig. 5 Maximized spanwise wavelength vs. Reynolds number 
for / from 2 kHz to 45 kHz for a Mach 2 flat plate boundary 
layer 

B.L. Length Scale (m) : 3.0E-09 - 2.4E-04 

U inf Velocity (m/s) : 590.065 - 590.065 

T inf Temperature (K) : 216.667 - 216.667 

\* */ 

% disturbance reference frequency in Hz % 

ranging from 153150 to 88310.4 for 1/d = 2 
or from 25524.9 to 14718.4 for 1/d = 12 

In addition to mean flow characteristic, a disturbance fre- 
quency range for unstable waves is also suggested. Based 
on this frequency range, an input file containing the fol- 
lowing may be used to find maximized growth rates with 
respect to the spanwise wavenumber i : 

// 

// Mach 2 flat plate 

// 

num_normal_pts =101 

mflow_f ilename = . /meanf low/mf low . m2fp" 

march ing_method_2d = along_re 

re = 500, final re = 7900, 

scep_re = 100, 

solut ion_type = local_eig_solution 


Title 

GasModel 

Reference Mach Number 
Prandtl Number 
No. of Stations 
X Coordinate Range (m) 
Re Range 


Mach 2 Flat Plate lod_max =25 


Perfect gas model 


1.99986 

0.7 

401 

0.000 - 
0.100 - 


f req_unit 
beta_uni t 

1.951 freq = 2e3, 
8000.00020e3, 25e3, 


= in_hertz_f req 
- wave_angle_beta 

4e3, 6e3, 8e3, 10e3, 12e3, 
30e3, 35e3, 40e3, 45e3, 


15e3, 
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beta = 13*60, 

qp_approx = true 

Figure 4 shows the resulting N-factor versus Reynolds 
number for various disturbance frequencies. The corre- 
sponding spanwise wavelength in millimeters is shown in 
Fig. 5. It indicates that the unstable spanwise wavelength 
is around 4 to 20 mm. The N-factor shown in Fig. 4 may be 
used for transition correlation or prediction directly. How- 
ever, it has been found that N-factor correlation based on 
individual modes rather than maximized modes represents 
the physics better and thus gives better transition correla- 
tion with respect to experimental data. Given the unstable 
frequency and spanwise wavelength range, further calcula- 
tions such as LST or linear PSE may be done for individual 
modes. Or if the initial disturbance amplitudes are known 
from receptivity calculations or experimental data, nonlin- 
ear PSE may be launched for a number of unstable waves 
to track transition onset. 



Fig. 6 Evolution of disturbance modal amplitudes for a Mach 
2 flat plate undergoing an oblique-mode breakdown 

As shown in Fig. 4, a disturbance with a frequency of 12 
kHz reaches an N-factor of 10 first. The associated max- 
imized spanwise wavelength is about 8 mm. We initiate 
a nonlinear PSE calculation dX R - 1500 with a pair of 
oblique modes with / = 12 kHZ and \ ~ mm (Fourier 
mode (1, 1) and (1, -1)). The initial amplitude is 0.1% 
for both modes. The resulting streamwise velocity ampli- 
tude versus Reynolds number is shown in Fig.6 for various 
Fourier modes. Mutual interaction of the oblique pair gives 
rise to the excitation and growth of the (0, 2) mode. Non- 
linear interaction of this triad then cascades energy into 
higher harmonics and eventually causes transition. The 
skin friction is plotted in Fig.7. The rapid rise near the 
end illustrates that the flow is transitional. More details 
concerning this oblique-mode breakdown may be found in 



Fig. 7 Skin friction coefficient vs. Reynolds number for a 
Mach 2 flat plate boundary layer undergoing an oblique-mode 
breakdown 

Chang & Malik.^ 

This example shows a typical work flow how LASTRAC 
may be used for stability calculations and transition predic- 
tions. For a given mean flow, the unstable frequency and 
spanwise wavelength range may be found with minimum 
effort and without too much prior knowledge of the flow. 
After identifying this range, N-factor correlations or pre- 
dictions using either LST or linear PSE may be performed. 
For a given initial amplitude, nonlinear PSE may also be 
used to track nonlinear development of the instability wave 
all the way until skin friction rise. 

Mach 6 Flared Cone 

Another test case is a Mach 6 flared cone for which ex- 
perimental measurements have been conducted in the Lan- 
gley Mach 6 2()-in tunnel.^’ The experimental model is a 

5 ’ straight cone for the first 10 in. followed by a flared 
section. At a free-stream Mach number of 6, the bound- 
ary layer edge Mach number inside the shock is about 
5.4. The dominant instability mode is hence second mode. 
The second-mode wave is most unstable when it is two- 
dimensional or axisymmetric. All second-mode calcula- 
tions presented here were done for axisymmetric waves. 
The asymmetric first-mode waves were computed by using 
different azimuthal wave numbers, ??, defined as the num- 
ber of waves along the azimuthal direction. 

Figure 8 shows quasi-parallel LST N-factors of various 
disturbance frequencies (ranging from 20 to 340 kHz) ver- 
sus the streamwise distance x for a unit Reynolds number 
of 2.89 X 10*' /ft. This case corresponds to the earlier Mach 

6 quiet wind-tunnel measurements with a unit Reynolds 
number of 2.81 x 10-' /ft. The N-factor value at the tran- 
sition onset is merely 3.8 under the conventional tunnel 
conditions, as compared to a value of about 7.8 under the 
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Fig. 8 Second-mode N-factors of various disturbance fre- 
quencies for the Mach 6 flared cone with a unit Reynolds 
number of 2.89 x lO” /ft at adiabatic wall conditions 



X (in) 

Fig. 9 Optimized first mode N-factor vs. Reynolds number 
for the Mach 6 flared cone, f = 60 kHz, also showing corre- 
sponding integer azimuthal wave number on the plot 

quiet tunnel conditions. 

Figure 9 depicts the N-factor and corresponding maxi- 
mized integer azimuthal wave number for a typical first- 
mode frequency of 60 kHz. The kink on the wave number 
curve near j' = 10 in. is due to the curvature discontinuity 
in the presence of the flare. The results show that the first 
mode is most unstable for n ranging from 10 to 50. Using 
this information, first-mode N-factors are calculated for a 
disturbance frequency ranging from 40 to 100 kHz and an 
azimuthal wave number from n = to n = 30 and the 
results are shown in Fig. 10. The most amplified azimuthal 


Fig. 10 First-mode N-factors of various disturbance frequen- 
cies and azimuthal wave numbers for a unit Reynolds number 
of 2.89 X 10^ /ft at an adiabatic wall condition 

wave number was found to be around 20 for most cases. 
The first-mode N-value at transition onset measured in the 
conventional (non-quiet) tunnel is only about 2 to 3 which 
is comparable to the second mode. 

Jet and Shear Layer 



y 


Fig. 11 Velocity and temperature profiles for a compressible 
shear layer 

LASTRAC may also be used for jet and shear layer insta- 
bility calculations. A shear layer with hyperbolic tangent 
velocity and temperature profiles, as shown in Fig. 1 1 , is 
used as a test case. The computed velocity and tempera- 
ture eigenfunctions are shown in Fig. 12. LASTRAC can 
be used to compute axi symmetric jets or vortex flows. The 
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Fig. 12 Velocity and temperature eigenfunction of an unsta- 
ble mode for a compressible shear layer at f = 8000 Hz 



Fig. 13 Growth rates vs. disturbance frequency for the Mach 
2.5 jet at X = 0.02 m for different azimuthal wave numbers 

test case is a Mach 2.5 jet. The computed growth rates 
versus disturbance frequencies are show in Fig. 13 for vari- 
ous azimuthal wave numbers ranging from 0 (axisymmetric 
mode) to 3. 

Thermal Control of Supersonic Crossflow Ins tability 

Recently, the concept of supersonic thermal laminar flow 
control has been proposed as an alternative to existing tech- 
nologies based on either natural laminar flow design or 
suction. Thermal control does not suffer from major draw- 
backs of shaping or suction and has the potential to offer an 
effective control with lower system penalties and fewer re- 
strictions on airframe design. This case is intended to study 


the effect of cooling on supersonic swept wing boundary 
layers. 

Viscosity is a strong function of temperature; it fol- 
lows that wall cooling directly affects the stability char- 
acteristics of viscous-mode instability waves. For gases, 
heating destabilizes while cooling stabilizes. Conversely, 
heating stabilizes a water boundary layer. Examples 
of viscous modes include the incompressible Tollmien- 
Schlicting (TS) wave and the first-mode wave in compress- 
ible boundary layers which is in fact a mixture of viscous 
and inviscid mode.^’ For inviscid modes, the effect of wall 
cooling is more intricate. The addition of a new general- 
ized inflection point due to wall cooling directly alters the 
stability characteristics. Mack"*- performed both inviscid 
and viscous stability calculations on first- and second-mode 
waves and concluded that in general, cooling stabilizes 
first-mode disturbances and it destabilizes higher modes 
such as a second-mode disturbance. 

For swept wing boundary layers, crossflow instability is 
the dominant instability mechanism. This inviscid insta- 
bility arises because of the presence of an inflection point 
in the crossflow profile in the direction perpendicular to the 
free-stream direction. As a result, for incompressible swept 
wing boundary layers, wall cooling has little effect on its 
growth. On the other hand, wall cooling may still affect su- 
personic crossflow instability in a significant way because 
of the presence of an additional generalized inflection point 
in the mean flow profiles. 



Fig. 14 Generalized inflection point for the streamwise ve- 
locity component for an adiabatic and cooled wall, Mach 2 
infinite swept wing boundary layer 

We use a Mach 2 flow over a 70® infinite swept wing 
boundary layer as the test example for thermal control. 
Before showing the stability results, we will examine the 
effect of cooling on the mean flow profiles first. The nec- 
essary condition for inviscid instability in compressible 
boundary layers is the presence of a generalized inflection 
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Adiabatic Wall 




Fig. 15 Generalized inflection point for the crossflow velocity 
component for an adiabatic and cooled wall, Mach 2 inflnite 
swept wing boundary layer 

point characterized by 


The generalized inflection quantity 
above boundary layer with an adiabatic and cooled wall 
{Tu:/Tadu = condition is shown in Fig.l4 and 15 

for the streamwise and crossflow profiles, respectively. As 
can be seen in Fig. 15, cooling simply moves the inflection 
point of the crossflow velocity profile toward the wall but 
does not remove it. In contrast (Fig. 14), two streamwise 
inflection points in the adiabatic wall case disappear as we 
apply cooling. It is interesting to note that the inflection 
points in the streamwise velocity profile are not present in 
the incompressible case. 

Our linear stability calculations show that similar to first- 
mode waves, cooling stabilizes both stationary and travel- 
ing low supersonic crossflow instability. The main differ- 
ence is that cooling is much more effective for first-mode 
waves that appear in 2-D boundary layers. For a Mach 2 
flat plate, Tu /Tadw = 0-95 is sufficient to reduce the N- 
factor significantly. However, in the present Mach 2 swept 
wing, a cooling of lujladw = 0.75 or lower is neces- 
sary to have a significant effect. Nevertheless, cooling still 
is an effective means to control supersonic crossflow in- 
stability. Figures 16 and 17 show the N-factor versus x/c 
for stationary crossflow instability with an adiabatic and 
partially cooled (with l\JTadu = 0.75 for .r/r < 0.1) 
wall, respectively. If we take A = 1 0 as the transition 
location, this partially cooled wall brings transition loca- 
tion from ,r/c = 0.52 to x/c = 0.41, a 15% reduction. 
The above results indicate that for supersonic crossflow 
instability, thermal control is a viable alternative to exist- 
ing technologies. More parametric studies are necessary to 


Fig. 16 N-factor vs. x/c for stationary crossflow for a Mach 
2, ^0" swept wing with an adiabatic wall condition 


TwTTadw » 0.75 up to x/c = .1 0 


Z 



Fig. 17 N-factor vs. x/c for stationary crossflow for a Mach 
2. 70 ^ swept wing with T„ ITad»r = 0 . 7-5 for x/c < 0.1 


draw a more concrete conclusion. 

Distributed Roughness Control 

This test case is concerned with the distributed rough- 
ness control proposed by Saric.^'^ The incompressible Ari- 
zona State University swept wing^^ boundary layer is first 
analyzed by using LST. Unstable stationary crossflow dis- 
turbances are identified by using the eigenvalue solution 
procedure. The unstable modes have a spanwise wave- 
length from 3 to 30 mm. Figure 18 shows some of the rep- 
resentative modes. In a distributed roughness control, we 
manually excite an early growing mode such as the 8-mm 
mode in this case using a distributed roughness near the 


13 

American Institute of Aeronautics and Astronautics Paper 





AIAA-2003-0974 


ASU Swept Wing 



Fig, 18 Important stationary crossflow modes for the incom- 
pressible Arizona State University swept wing experiment 


Velocity Perturbation Amplitude 



Fig. 19 Nonlinear evolution of the 12-mm and 24-mm modes 
with and without the presence of the controlled 8-mm mode 
for the incompressible ASU swept wing experiment 

leading edge. This early growing mode will grow to a large 
amplitude near the leading edge and its growth modifies the 
mean flow to an extent that the normal growth of the most 
amplified mode with a wavelength of 12 mm would be sup- 
pressed. In the mean time, the 8-mm mode rapidly enters 
the decaying stage downstream because of the nature of the 
short wavelength modes. Eventually, all unstable modes 
involved die out and laminar flow is maintained before the 
end of the chord. Nonlinear PSE was run for two cases. 
The uncontrolled case was run with an initial amplitude of 
about 0.01 % for the 12-mm mode. For the controlled case, 
we introduce an additional 8-mm mode with an initial am- 


plitude of about I . As shown in Fig. 1 9, the uncontrolled 
case reaches an amplitude of about 1 -VX near the end of 
simulation. On the other hand, the controlled case has a 
final amplitude less than 2% for both 8-mm and 12-mm 
modes (as well as its harmonic 24-mm mode). Transition 
in this configuration is related to secondary instability of 
the highly nonlinear stationary crossflow modes. Smaller 
amplitude in the controlled case would prevent secondary 
instability and thus transition from happening. These cases 
were run with 24 Fourier modes (M — 0, A' = 24). 


Mach 2, 70 degree Swept Wing 



Fig. 20 Important stationary crossflow modes for the Mach 
2, 70"^ swept wing boundary layer 



Fig. 21 Nonlinear evolution of the 12-mm and 24-mm modes 
with and without the presence of the controlled 8-mm mode 
for the Mach 2, 70'^ swept wing boundary layer 

To investigate whether a distributed roughness control 
would work for supersonic boundary layers, we study the 
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Veloc^ Perturbation Amplitude 



Fig. 22 Nonlinear evolution of the 12-nim and 24-nun modes 
with and without the presence of the controlled 6-mm mode 
for the Mach 2, 70 swept wing boundary layer 

Mach 2, 70'^ swept wing configuration used previously for 
thermal control again. Figure 20 shows the N-factor for the 
most amplified along with other representative modes. The 
most amplified wave is at — 12 mm and both 6-mm and 
8-mm waves appear to be a good candidate for the control. 
Note that the 24- mm mode also has a large N-factor at the 
end. For the uncontrolled calculation, we excite both the 
12-mm and 8-mm modes with an initial amplitude of 1 - 

6. The fundamental spanwise wavelength is 24 mm and 32 
modes (A' = .32) are kept in the Fourier series. 

The first control is to utilize the 8-mm mode. We in- 
crease the initial amplitude of the 8-mm mode an order 
of magnitude larger to o.r - 5. The results are shown in 
Fig. 21. For the uncontrolled case, the 8-mm mode grew 
early on but it was overtaken by the 1 2- mm mode at around 
xjc = 0.22, the harmonic 24-mm mode also grew to a 
quite large amplitude near the end of the calculation. For 
the controlled case, the 8-mm mode was dominating early 
on then started to decay not too far downstream. As a result, 
both 12-mm and 24-mm modes were suppressed. The final 
amplitude for the controlled case is at least one and half 
orders of magnitude smaller than the uncontrolled case. 

For the second control, we introduce the 6-mm mode as 
the (0, 4) mode in the nonlinear calculation. The initial am- 
plitudes of the 12-mm and 8-mm modes are set to 1 .f - 6 
and the control mode 6-mm has an amplitude of 5.f - 5- 
Figure 22 shows the results. The 6-mm control mode grew 
early on and was dominating for xjc < 0.18, It then de- 
cayed rapidly. Similar to the 8-mm control case discussed 
previously, the growth of 12-mm and 8-mm modes was 
suppressed by the presence of the large amplitude 6-mm 
mode. However, the rapid growth of the 247???7/ mode at 
the end appears to be unavoidable. This 24-mm mode ex- 
cited by nonlinear interaction of the 6-mm and 8-mm ((0, 4) 


and (0, 3) modes) as well as 8-mm and 12-mm ((0, 3) and 
(0, 2) modes) grew very rapidly from the beginning. Fur- 
thei-more, for xjc > 0.2.3 where other modes are decaying, 
this 24-mm mode begins to dominate due to its large linear 
growth rate (see Fig.20). Therefore, the 6-mm control is not 
as effective as the 8-mm one as shown in Fig. 21. Depend- 
ing on the initial amplitude, in this case, the introduction 
of the 6-mm mode may or may not achieve laminar flow 
control because of the presence of the ver>^ unstable longer 
wavelength mode. 

This example shows that the selection of the control 
mode is crucial to the success of distributed roughness 
control. For the latter supersonic swept wing case, lin- 
eally, the unstable stationary crossflow modes cover a large 
spanwise wave number range, which is a characteristic of 
large Reynolds number configurations. As a result, non- 
linear interaction of the control mode and the most am- 
plified mode may give rise to a longer wavelength mode 
that would eventually dominate downstream. When this 
happens, using distributed roughness with a single span- 
wise wavelength would fail according to our nonlinear PSE 
results shown here. However, introducing multiple-scale 
distributed roughness may still work. Further parametric 
study is necessary to draw a more concrete conclusion. 
The results presented here also demonstrate the importance 
of parametric study using nonlinear PSE for a nonlinear 
control concept such as the distributed roughness control. 
LASTRAC provides one such tool for future design stud- 
ies. 

Incompressible Subharmonic Breakdown 



Fig. 23 Nonlinear evolution of selected Fourier modes for 
subharmonic breakdown to turbulence in a Blasius boundary 
layer 

The last test case is the well-known incompressible sub- 
harmonic breakdown phenomena. The fundamental mode 
is a 2-D mode (mode (2.0)) with a frequency of f = 
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Fig. 24 Skin friction coefficient versus Reynolds number for 
the subharmonic breakdown in a Blasius boundary layer 

1.24 X The subharmonic mode (mode (1, 1)) has 

a nondimensional (normalized with the boundary layer 
length scale) spanwise wavenumber of about /i = 0.14. 
Initial amplitudes are about OJM and 0.005% for the fun- 
damental and subharmonic modes, respectively. We use 
six Fourier modes in both temporal and spanwise direc- 
tions (M — 6, A = 6). Figure 23 shows the u velocity 
amplitude evolution for the fundamental, subharmonic and 
some representative modes. Secondary instability takes 
place early on; and nonlinear interaction results in exci- 
tation of many harmonics. Near the end of the simulation, 
the fundamental mode is forced to grow again (back reac- 
tion) as the subharmonic mode saturates. The skin friction 
for the perturbed nonlinear PSE solution is plotted against 
the Blasius (laminar) solution in Fig. 24. A rapid rise of 
skin friction is evident. 

Summary 

Predicting transition onset remains a daunting task even 
with the state-of-the-art computing facilities. At the lowest 
fidelity, transition may be predicted by using a prescribed 
N-factor and a simple linear stability or parabolized sta- 
bility theory. Depending on the accuracy requirement, the 
N-value may come from empirical correlations with exist- 
ing wind tunnel or flight experiments, or from a commonly 
used value such as 10. The major problem of this approach 
lies in the N-value itself. According to past experiences, 
the transition N-value may vary from as small as 2 or 3 
for a noisy facility to as high as 15 or 20 in flight. Thus, 
for a new configuration, determining the N value is itself a 
difficult task. Despite this uncertainty, N-factor correlation 
remains the most viable method for transition prediction 
due to its simplicity. LASTRAC provide both LST- and 
PSE-based N-factor correlation capability. 

For a given mean flow, a possible unstable frequency 


range is suggested by LASTRAC and by using the max- 
imizing N-factor option, users may obtain an unstable 
parameter range of disturbance frequency and spanwise 
wavelength with minimum effort and little prior knowledge 
of the mean flow under consideration. Further linear calcu- 
lations using either quasi-parallel LST or nonparallel PSE 
may be launched to obtain an envelop of N-factors formed 
by a broad range of unstable modes. The N-factor envelope 
then can be used for transition correlations or predictions. 

LASTRAC also provides the capability to compute dis- 
turbance evolution based on an absolute amplitude. Non- 
linear PSE calculations may be performed for a number 
of unstable modes with a finite amplitude. We show sev- 
eral test cases in which the boundary layer is perturbed 
with a given amplitude and transition is captured by us- 
ing the nonlinear PSE option. Coupled with the receptivity 
model which will be incorporated in the near future, LAS- 
TRAC offers a transition prediction tool that may be used to 
compute transition onset without any modeling or N-factor 
assumptions. 

It is also demonstrated that LASTRAC may be used for 
parametric studies of several supersonic laminar flow con- 
trol concepts. We presented two such techniques, thermal 
and distributed roughness control. The test cases shown in 
this paper cover a broad range of flow configurations. In ad- 
dition to the traditional N-factor method, LASTRAC offers 
a comprehensive set of options based on the state-of-the-art 
numerical methods that may be used for the stability calcu- 
lations and transition predictions in an integrated fashion. 
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